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ABSTRACT 

We compare deep VLA imaging of the total intensity and linear polarization of the 
inner jets in the nearby, low-luminosity radio galaxy 3C 31 with models of the jets as 
intrinsically symmetrical, decelerating relativistic flows. We show that the principal 
differences in appearance of the main and counter-jets within 30arcsec of the nu- 
cleus can result entirely from the effects of relativistic aberration in two symmetrical, 
antiparallel, axisymmetric, time-stationary relativistic flows. We develop empirical pa- 
rameterized models of the jet geometry and the three-dimensional distributions of the 
velocity, emissivity and magnetic-field structure. We calculate the synchrotron emis- 
sion by integration through the models, accounting rigorously for relativistic effects 
and the anisotropy of emission in the rest frame. The model parameters are optimized 
by fitting to our 8.4-GHz VLA observations at resolutions of 0.25 and 0.75arcsec 
FWHM, and the final quality of the fit is extremely good. The novel features of our 
analysis are that we model the two-dimensional brightness distributions at large num- 
ber of independent data points rather than using one-dimensional profiles, we allow 
transverse as well as longitudinal variations of velocity, field and emissivity and we 
simultaneously fit total intensity and linear polarization. 

We conclude that the jets are at «52° to the line of sight, that they decelerate and 
that they have transverse velocity gradients. Their magnetic field configuration has 
primarily toroidal and longitudinal components. The jets may be divided into three 
distinct parts, based not only on the geometry of their outer isophotes, but also on 
their kinematics and emissivity distributions: a well-collimated inner region; a flaring 
region of rapid expansion followed by recollimation and a conical outer region. The 
inner region is poorly resolved, but is best modelled as the sum of fast (0.8 - 0.9c) and 
much slower components. The transition between inner and flaring regions marks a 
discontinuity in the flow where the emissivity increases suddenly. The on-axis velocity 
stays fairly constant at «0.8c until the end of the flaring region, where it drops abruptly 
to «0.55c, thereafter falling more slowly to «0.25c at the end of the modelled region. 
Throughout the flaring and outer regions, the velocity at the edge of the jet is ~0.7 of 
its on-axis value. The magnetic field in the flaring region is complex, with an essentially 
isotropic structure at the edge of the jet, but a more ordered toroidal+longitudinal 
configuration on-axis. In the outer region, the radial field vanishes and the toroidal 
component becomes dominant. We show that the emissivity and field structures are 
inconsistent with simple adiabatic models in the inner and flaring regions. We suggest 
that the discontinuity between the inner and flaring regions could be associated with 
a stationary shock structure and that the inferred transverse velocity profiles and field 
structure in the flaring region support the idea that the jets decelerate by entraining 
the external medium. We demonstrate the appearance of our model at other angles 
to the line of sight and argue that other low-luminosity radio galaxies resemble 3C 31 
seen at different orientations. 

Key words: galaxies: jets - radio continuum:galaxies - magnetic fields - polarization 
- MHD 
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1 INTRODUCTION 

The flow parameters of jets in extragalactic radio sources 
have hitherto proven difficult to determine because of the 
absence of unambiguous diagnostics. Most progress has been 
made in the estimation of velocities, particularly where these 
are thought to be r elativistic. The idea th at jets in low- 
luminosity, i.e. FRI (Fanaroff & Riley 1974), radio galaxies 



have relativistic speeds rests on five main arguments: 



(i) evidence for relativistic motion on parsec scales in BL 
Lac objects, coupled with the hypothesis that they are FRI 
r adio galaxies observed at small angles to the line of sight 
( |Urry fc Padovani 1995| ); 

(ii) measurement of a pparent proper motions with s peeds 
up to at least c in M87 (Biretta, Zhou & Owen 1995); 



(iii) modelling of relativistic flows, which demonstrates 
the feasibility of decelera tion on kiloparsec scales in realis- 



tic galactic atmospheres (Bickncll 1994 ; 
Bowman, Leahy & Komissarov 1996); 



Komissarov 1994 



(iv) the interpretation of correlated depo larization asym- 



metr y and jet sidedness in FRI sources (Morganti et al 



1997) as a conseq uence of Dop pler beaming and foreground 



Faraday rotation (Laing 1988); 

(v) observations of brightness and width asymmetries in 
FR I jets which decrease with distance from the nucleus and 
are correlated with fractional core flux, implying tha t they 



decelerate and are faster on-axis than at their edges ( Laing 



1992 



1996; Hardcastle et al. 1997; Laing et al. 199S) 



This paper reports a detailed study of the applicability 
of decelerating relativistic jet models to the well-resolved 
kiloparsec-scale structures in the nearby FRI radio galaxy 
3C31. Our intention is to deconvolve the emission mecha- 
nism (synchrotron radiation from a magnetized relativistic 
flow) from the radio data, without embodying specific pre- 
conceptions about the poorly known internal physics. We 
construct sophisticated three-dimensional models of the ef- 
fects of relativistic aberration on the appearance of intrinsi- 
cally symmetric magnetized jets and we fit these models to 
the observed total and polarized intensity distributions in 
3C31. The aim is to derive robust estimates of the velocity 
field, the emissivity (combining relativistic particle density 
and magnetic field strength) and the three-dimensional or- 
dering of the magnetic field (purely geometrical factors in- 
dependent of its strength). We regard this as a necessary 
first step: realistic physical models capable of being com- 
pared with observations are not yet available. We are able 
to reproduce many of the observed features of the jets and 
we conclude that our models now provide a way to obtain 
key constraints on the intrinsic properties of extragalactic 
radio jets. 

Section ^ describes new VLA imaging of the jets in 
3C 31 that provides a high-quality data set suitable for de- 
tailed fitting by our models. Section ^| first reviews the prin- 
ciples underlying the relativistic jet models, and then goes 
on to describe how we adjust their parameters to fit the ra- 
dio intensity and polarization data. Section ^ critically dis- 
cusses the model fits and reviews the main features of the 
inferred jet velocity field, magnetic structure and emissivity 
distribution. Section ^ discusses more general implications, 
including: specific problems associated with reproducing the 
properties of the jets in the region closest to the galactic nu- 



cleus; the reasons for the sudden onset of deceleration; ev- 
idence for interaction between the jet and the surrounding 
medium and the applicability of adiabatic models. 

Section ^ illustrates how the jets in 3C 31 should appear 
if orientated at other angles to the line of sight and outlines 
the applicability of the model to other FRI sources, includ- 
ing those orientated at a small angle to the line of sight. 
Section [?] summarizes our conclusions regarding the kine- 
matics, emissivity and field structure of the three distinct 
regions of the jet and their implications for future work. 

Throughout this paper, we adopt a Hubble constant Ho 
= 70 km s _1 Mpc _1 . 



2 OBSERVATIONS AND IMAGES 

2.1 Choice of source: the assumption of intrinsic 
symmetry 

Our key assumption is that the bases of the two jets are 
intrinsically identical, antiparallel, axisymmetric stationary 
flows, and that the apparent differences between them result 
entirely from relativistic aberration. 

This is an approximation in two important respects: we 
ignore small-scale structure in the jets, and we assume that 
any intrinsic or environmental asymmetries (clearly domi- 
nant on the largest scales in many objects) are small com- 
pared with relativistic effects close to the nucleus. For any 
individual source, it is difficult to be sure that this is the 
case, particularly if the asymmetry persists on all scales 
Ther e are also a few sources (e.g. 0755+379 



Bondi et al 



200C) in which the counter-jet appears much wider than the 
main jet, an effect which cannot be produced by relativis- 
tic beami ng. We argue that these cases are rare (Laing et 



al. 1999) and easy to recognize. If the jet/counter-jet ra- 



tio decreases with distance from the nucleus, approaching 
unity on large scales, then relativistic effects probably dom- 
inate, as most plausible intrinsic or environmental mecha- 
nisms would generate asymmetries which stay roughly con- 
stant or even increase with distance. A statistical study of 
a complete sa mple of FRI sourc es with jets selected from 
the B2 sample ( fr^aing et al. 199*9 ) suggests that the median 
asymptotic jet/counter-jet ratio on large scales ~1.1 for the 
weaker sources, so the assumption of intrinsic symmetry for 
the jets is generally reasonable. We ensure that this assump- 
tion is self-consistent by choosing an object whose jets are 
straight, with similar outer isophotes on both sides of the 
nucleus. 

We also require that the jets are bright, allowing imag- 
ing with high signal-to-noise ratio in total intensity and lin- 
ear polarization, and that any effects of Faraday rotation can 
be accurately corrected. This led us to the choice of 3C 31 
as the first source to model. 



2.2 3C31 

3C 31 is an F R I radio galaxy that has long been known (e.g. 
Burch 1977) to have a resolved bright jet and counter-jet 



with a strong initial brightness asymmetry, also present on 
parsec scales ( Lara et al. 1997 ) , a hi gh degree of linear po - 
larization at centimetre wavelengths ( Fomalont_et_aL_r98C|) 
and a significant depolarization asymmetry (Burch 1979j; 
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Strom et al. 1983 ). It is identified wit h the dusty elliptical 
galaxy NGC 383 ( Martel et al. 1999[ ) that is the bright- 
est me mber of a rich group of galaxies (Arp 1966; Zwicky 
Kow^l 196£) with an exte nsive hot intra-group medium 
BShringer 1999[ ). 



(Komossa . 



We adopt a redshift of 0.0169 for the galaxy (the mean 



of v alue s f rom Smith et al. 2000 , Huchra, Vogeley fc Geller 
1999] aJid IDe Vaucouleurs et al. 1991), so the projected 



linear scale of the images is 0.34 kpc/arcsec. 



Table 1. Journal of observations 



Configuration 



A 
B 
C 
D 



Frequency 
MHz 

8460 
8440 
8440 
8440 



Date 



1996 Nov 12 
1994 Jun 6, 14 

1994 Dec 4 

1995 Apr 28 



Integration 
time (min) 

606 
818 
242 
69 



2.3 VLA observations and data reduction 



The observations (Table j^j were made using all four configu- 
rations of the VLA at centre frequencies of 8.46 or 8.44 GHz 
and a bandwidth of 100 MHz (the slight difference in fre- 
quencies between configurations has no measurable effect). 
They were reduced in the AIPS software package using stan- 
dard self-calibration and imaging methods, with one major 
exception, viz. an iterative technique used to combine data 
from the different VLA configurations when the compact 
core had varied significantly between observations. In order 
to make the best possible images, we needed to adjust for 
these changes in the compact core, as well as for slight incon- 
sistencies in the amplitude calibration of the four observing 
runs. 

We therefore adopted the following procedure: 

(i) Image and clean the data from the widest (A) con- 
figuration, ensuring that the compact core is centred on a 
map piper: 



(ii) dftlf-calibiate, initially adjusting uiihy tha phases; then 

the amplitudes, until the best image is obtained. 

(iii) Use the clean components from this image as a 
model for phase-only self-calibration of the next widest (B) 
configuration. 

(iv) Image both datasets at the same resolution in order 
to measure the flux density of the core. 

(v) Adjust the (u,v) data for the larger configuration by 
adding or subtracting the appropriate point component to 
equalize the core flux densities. 

(vi) Concatenate the two (u,v) datasets, image, CLEAN 
and self-calibrate phases and amplitudes as in steps (i) and 
(ii). 

(vii) Split the datasets apart again and check that the 
core flux densities are consistent. If not, repeat steps (v) 
and (vi). 

(viii) Add further VLA configurations using steps (iii) to 
(vii). 

Two sets of images in Stokes /, Q and U were made from 
the combined four-configuration data set, one with full reso- 
lution (Gaussian FWHM 0.25 arcsec) and the other tapered 
to give a FWHM of 0.75 arcsec (Table ||). Both the Maxi- 
mum Entropy and CLEAN algorithms were used to compute 
deconvolved I images. The compact core was subtracted 
from the data before Maximum Entropy deconvolution, and 
added in again afterwards. All images were restored with 
the same truncated Gaussian beam. 

The result of differencing the Maximum Entropy and 
CLEAN / images was a high-frequency, quasi-sinusoidal rip- 
ple of near zero mean whose amplitude increased with sur- 
face brightness. This artefact clearly originated in the CLEAN 



Table 2. Images and rms noise levels. 



Resolution rms noise level 
(arcsec) (mJv / beam area) 

I Q/U 



0.25 
0.75 



5.5 
6.9 



6.1 
5.5 



image and such ripples are indeed known t o be characteris - 
tic of instabilities in the CLEAN algorithm (Cornwell 1983). 
There was no evidence for any differences between the two 
images on larger scales. We therefore use only the Maximum- 
Entropy / images in what follows. The Q and U i mages were 
CLEANed. A first-order correction for Ricean bias (Wardle & 



Kronberg 1974) was made when deriving images of polar- 



ized intensity. Apparent magnetic field directions were de- 
rived from images of Q and U corrected for Faraday rotation 
using results from a six-frequency analysis of the polarime- 
try of this region at a resolution of 1.5 arcsec FWHM to be 
published elsewhere. The maximum correction is «9° and 
Faraday depolarization is negligible. 

The resulting images are almost noise- limited (Table ^), 
and the excellent (u,v) coverage and signal-to- noise ratio 
allow a good representation of jet structures on a wide range 
of scales. 



2.4 Source description and selection of regions to 
be modelled 

The left panel in Fig. [j] shows the large scale structure of 
3C31 derived from lower-resolution observations at 1.4 GHz 
that will be published elsewhere. The jet and counter-jet 
bend on large scales and form two extensive sinuous plumes. 
The inner part of our full-resolution 8.4-GHz image is shown 
in the right panel. Both jets are well collimated and dim 
for the first 2.5 arcsec from the nucleus, flare and brighten 
between 2.5 and 8 arcsec and th en recollimate, as noted from 
earlier VLA imaging at 5 GHz ( [Fomalont et al. 1980| ). Our 



assumption of intrinsic symmetry requires that we restrict 
our analysis to the straight regions of the jets shown in the 
rectangular box in the right panel of Fig. |l|. Within this 
area, the outer isophotes in the main and counter-jets are 
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Figure 1. Montage showing the large-scale structure and jets of 3C31. Left panel: VLA 1.4-GHz image of a 15 arcmin (300 kpc) North- 
South field at 5.5arcsec (1.9 kpc) resolution. Right panel: VLA 8.4-GHz image of an approximately 2 arcmin (40 kpc) North-South field 
at 0.25 arcsec (85 pc) resolution. The rectangle within the right panel shows the relatively straight segment of the jets that we have 
chosen to model. 

5 10 15 

indeed very similar. The emission from this region is well- 
resolved and bright enough to provide strong constraints on 
the velocities both along and transverse to the jet at about 
1300 independent locations. 

The systematic differences in brightness and polariza- 
tion structure between the jets provide essential clues to the 
orientation and velocity field. An image of sidedness ratio, 
constructed by dividing the I image by a copy of itself ro- 
tated through 180° (in the sense main jet/counter-jet) is 
shown in Fig. ^ The ratio is ~ 5 close to the nucleus and 
has a mean of ~ 13 (with erratic fluctuations) from 2.5 to 
6 arcsec. It drops rapidly between 6 and 8.5 arcsec, there- 
after falling smoothly to ~ 1 at the end of the modelled 
region. The ratio at the jet edge is almost always lower than 
the on-axis value, but significantly exceeds unity except at 
distances from the nucleus > 20 arcsec. 




Figure 2. The observed jet/counter-jet brightness ratio (sided- 
ness) at a resolution of 0.75 arcsec, from the 8.4-GHz observations. 
This was constructed by dividing the / image by a copy of itself 
rotated through 180° and is in the sense main jet/counter-jet. 



Fig. g shows the degree of polarization and inferred 
magnetic field orientation over the inner ±27 arcsec of both 
jets at 0.75 arcsec resolution. The predominant pattern of 
the apparent magnetic field directions is to be perpendicular 
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Figure 3. Vector plot showing the degree of polarization, p, and apparent ma gnetic field direction superposed on a total intensity grey 
scale at a resolution of 0.75 arcsec. The three jet regions defined in Section |3.2| are indicated by vertical bars. 



to the jet axis near the centre lines of both jets, and parallel 
to the outer isophotes (with high degrees of polarization) at 
both jet edges. A notable exception to this trend occurs be- 
tween 6 and 10 arcsec from the nucleus at the edges of both 
jets, where the degree of polarization is very small. There is 
a polarization asymmetry on small scales: in the bright base 
of the main jet, the apparent magnetic field lies alo ng the jet 
axis (as is usual in FRI jets; Bridle & Perley 1984), whereas 



in the equivalent part of the counter-jet, it is transverse. In 
the outer regions, the degree of transverse-field polarization 
on-axis is significantly higher in the counter jet. 

The jets clearly develop some internal structure (arcs 
and non-axisymmetric knots) before they reach the large- 
scale bends (Fig. hi) . Similar structures are seen in linear 
polarization (Fig. H): most prominently, the apparent mag- 
netic field lies parallel to the intensity ridge lines in the arc- 
like structures seen on the total intensity images, produc- 
ing a vortex-like appearance in the magnetic vectors. These 
features, which account for some tens of per cent of the 
brightness locally, cannot be modelled assuming that the 
flow is smooth, axisymmetric and stationary. Together with 
the large-scale bends in the jet, they effectively set the limits 
to which we can fit the data. 



3 THE MODEL 
3.1 Assumptions 

Our key assumption is that the bases of the two jets are in- 
trinsically identical, antiparallel, axisymmetric, stationary 



(this tu rns out to be the case at one special location: see 
Section 3.3). 

We allow both for lon gitudinal decelera tion as inferred 
from the B2 sample data (Laing et al. 1999) and for trans- 
verse velocity structure. The latter is required for two rea- 
sons. First, one class of model for relativistic jet deceleration 
invokes entrainment of the interstellar medium from the host 
galaxy across a boundary layer. A transverse velocity vari- 
ation allows for the possibility that the relativistic particles 
near the edges of the jet move down the jet more slowly than 
those on the jet axis. Second, because the outer isophotes of 
the jets in FRI sources are usually more sym metric across 
the nucleus than those close to the jet axis ( Laing 199^ , 



199C ; Hardcastle et al. 1997 



Laing et al. 199S ) , a transverse 



flows. We model the jets using simple parameterized ex- 



velocity variation is generally required to fit well-resolved 
jet brightness distributions. In order to quantify this effect, 
we consider t wo possible transverse structures. In the first 
(Laing 1993), a central fast spine with no transverse varia- 
tion of velocity or emissivity is surrounded by a slower shear 
layer with gradients in both variables. In the second case, 
there is no distinct spine component, and the jet consists 
entirely of a shear layer with a truncated Gaussian trans- 
verse variation in velocity. We will refer to the two types as 
spine/shear-layer (SSL) and Gaussian models, respectively. 

The significant linear polarization observed requires 
an anisotropic magnetic field. We assume that it is disor- 
dered on small scales, with negligible mean, and that the 
anisotropy is introduced by shear and compression. We con- 
sider large-scale ordering of the magnetic fields to be un- 
likely. The simplest ordered fields capable of generating the 
observed polarization (single helices) produce large c hanges 
in em i ssion across the jets, which are not observed ( Laing 
1981; Clarke, Norman fc Burns 1989] ). More complex or- 



pressio ns for the variables which determine the synchrotron dered configuratio ns, such as those proposed by ^<onigl & 



emissio l - velocity fields, emissivity variations and intrinsic Choudhuri (1985), cannot be as easily dismissed on obser- 



magnetic field structures - and determine the free parame- 
ters of these expressions by fitting to the observed images. 

We assume that the flow is laminar and that there are 
no discontinuous changes of direction. If there is a turbu- 
lent velocity component (as theoretical models suggest; Sec- 
tion |5.2| ), then our technique will determine an average bulk 
flow speed, weighted by the distribution of Doppler beaming 
factors in a given region. We also assume that the variations 
of velocity and emissivity are continuous and smooth un- 
less our fitting procedure explicitly requires discontinuities 



vational grounds, but the presence of a significant ordered 
longitudi nal component is ruled out by flux co nservation ar- 
guments (Begelman, Blandford & Rees 1984). In any case, 



our conclusions on the relative magnitudes of the field com- 
ponents would not be seriously affected (indeed, our calcu- 
lations would be unchanged if one of the three field compo- 
nents is vector-ordered). We quantify the anisotropy using 
the ratios of the rms field components along three orthogo- 
nal directions. 

The spectrum of the jets between 1.4 and 8.4 GHz at 
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a resolution of 1.5arcsec FWHM is accurately described by 
a power law with a spectral index a = 0.55 (S v oc v~ a ~). 
The emission is therefore taken to be optically thin (we 
do not attempt to model the partially self-absorbed parsec- 
scale core). The corresponding electron energy distribution 
is n(E)dE = n E' {2a+1) dE = n E~ 21 dE. We assume an 
isotropic pitch-angle distribution relative to the field, so the 
degree of polarization, p = (U 2 + Q 2 ) 1 ^ 2 /I , has a maximum 
value of po = (3a + 3)/(3a + 5) = 0.70. 



3.2 Geometry 

We define to be the angle between the jet axis and the 
line of sight, z is a coordinate along the jet axis with its 
origin at the nucleus, x is measured perpendicular to the 
axis, r = (x 2 + z 2 ) 1 ' 2 is the distance from the nucleus and 
(f> is an angle measured from the jet axis (x = ztan</>). The 
first step in our procedure is to define functional forms for 
the outer surfaces of the jets and for the flow streamlines. 
The latter inevitably involves some guesswork, to be justified 
post hoc by the quality of the model fit. Inspection of the 
outer isophotes shows that the jets can be divided into three 
regions: 

(i) Inner (0 - 2.5 arcsec): a cone, centred on the nucleus, 
with a half-opening angle of 8.5 degrees. 

(ii) Flaring (2.5 - 8.3 arcsec): a region in which the jet 
initially expands much more rapidly and then recollimates. 

(iii) Outer (8.3 - 28.3 arcsec): a second region of coni- 
cal expansion, also centred on the nucleus, but with a half- 
opening angle of 16.75 degrees. 

All dimensions given above are as observed, i.e. projected on 
the plane of the sky. This pattern of an initially narrow base 
and a rapid expansion followed by recollimation is general 
in FRI jets ( Bridle fc Perley 1984 ). In what follows, we use 
subscripts i, f and o to refer to quantities associated with 
the inner, flaring and outer regions. We refer to the inner 
and outer boundaries separating the regions by subscripts 
1 and 0. The inner b ound ary is the flaring po int defined by 
Parma et al. (1987) and Laing et al. (1999), and we also 



Table 3. Definitions of streamline indices for inner and outer 
regions. 



use this term. 

Guided by the shape of the outer isophotes, we assume 
that the flow in the inner and outer regions is along straight 
lines passing through the nucleus. Our general approach is 
to devise simple analytical functions to describe the flow 
in these regions, and then to interpolate across the more 
complex geometry of the flaring region in such a way as to 
preserve continuity. Families of streamlines are parameter- 
ized by the streamline index s, which varies from at the 
inner edge of a component (spine or shear layer) to 1 at the 
outside edge. In the inner and outer regions, the streamlines 
make constant angles <j>i and (j> with the jet axis. We define 
£i and £ to be the half-opening angles of the jet in the inner 
and outer regions, and £i, £<, to be the corresponding angles 
for the spine, s is defined in terms of these angles in Table [| 

We require continuity of the streamlines and their first 
derivatives with respect to z across the flaring region. The 
simplest functional form that satisfies these constraints and 
fits the outer isophote shape for s — 1 in the shear layer is: 

x = ao(s) + ai(s)z + a2(s)z 2 + as(s)z 3 



Model 



Inner 
(conical) 



Outer 
(conical) 



SSL spine 
SSL shear layer 
Gaussian 



Ci s <Po = CoS 

Ci + (6 - CO* <f>o = Co + (Co - Co)s 



For each streamline, the values of ao(s) - 03(5) are deter- 
mined uniquely and in analytic form by the continuity con- 
ditions. The natural boundaries between regions are then 
spherical, centred on the nucleus at distances n and ro 
and therefore perpendicular to the streamlines. Fig. ^ shows 
sketches of the assumed geometry for the SSL model (the 
equivalent for the Gaussian model is essential identical, but 
with the spine removed). 

In order to describe variations along a streamline, we 
use a coordinate p, defined as: 



p — r (inner region) 
p = ri + (r — n)- 



- n cos <f>i(s) 



p = 



ro cos (j> (s) — ri cos 4>i(s) 
(flaring region) 
r (outer region) 



p is monotonic along any streamline and varies smoothly 
from n to ro through the flaring region (p — r — z on the 
axis). This allows us to match on to simple functional forms 
which depend only on r. 

The functions defining the edge of the jet are con- 
strained to match the observed outer isophotes and are fixed 
in a coordinate system projected on the sky. Their values in 
the jet coordinate system then depend only on the angle to 
the line of sight. The outer edge of the spine in SSL models 
is not constrained in this way, and the relevant parameters 
may be varied in order to obtain a good fit to the data. 

In what follows we will refer to streamline coordinates 
defined by longitudinal (along a streamline), radial (out- 
wards from the axis) and toroidal orthonormal vectors. 



3.3 Parameter variations 

Our approach to parameterizing the variations of velocity, 
emissivity and field ordering along the jets is to specify val- 
ues at four standard locations: inner jet, inner boundary 
(just inside the flaring region) , outer boundary and an arbi- 
trary fiducial point in the outer region. There is insufficient 
information to constrain any variations along the inner re- 
gion, so constant values are assumed there. We allow discon- 
tinuities in most variables at the inner boundary, as there 
is unambiguous evidence for an abrupt change in the emis- 
sivity, at least, at this position. By contrast, all quantities 
vary continuously through the flaring and outer regions. We 
specify their values at the three locations, together with any 
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Figure 4. Geometry of the spine/shear-layer model, showing the 
inner, flaring and outer regions in the plane containing the jet 
axis. The thick full curves represent the edge of the jet, the bound- 
aries between regions are represented by thin full curves and the 
s = 0.5 streamlines for the spine and shear layer are drawn as 
dashed curves, (a) The entire modelled region; (b) the base of the 
jet on a larger scale, showing the boundary surfaces at distances 
of n and ro from the nucleus. The Gaussian model is essentially 
the same, but with the spine component removed. 



uniformly. We chose simple functional forms for [3 p (p) to 
satisfy these requirements (see Table ^). The constants 60 
- &2, Co and ci are chosen to match specified velocities at 
the inner and outer boundaries, and at an arbitrary fiducial 
point in the outer region. In addition, we require continu- 
ity of velocity and acceleration across the outer boundary 
so the constants are uniquely determined. These continuity 
conditions are not strictly necessary for the calculation de- 
scribed here, but are physically reasonable and essential for 
the adiabatic models that we discuss elsewhere. 

In the flaring and outer region, /3 S = 1 in the spine 
for SSL models, dropping linearly with s from 1 at the 
spine/shear layer interface to a minimum value at the edge 
of the jet. For Gaussian models, f3 s is a truncated Gaussian 
function. The fact that the sidedness ratio at the edge of the 
jets exceeds unity over most of the modelled area (Fig. ^|) 
means that the fractional velocity at the edge is significantly 
greater than zero. In both classes of model, this minimum 
fractional velocity, v(p), is allowed to vary along the jet (Ta- 
ble |). 

In the inner region, we found that we could not obtain 
satisfactory fits with linear or Gaussian transverse velocity 
profiles. The data required a mixture of fast and slow mate- 
rial, without much at intermediate velocity. Given the poor 
transverse resolution, we took the simple approach of as- 
signing a single fractional velocity f3 s = V\ to material in 
the "shear layer" (there is actually no shear). This means 
that the velocities of the spine and shear layer in the SSL 
models are decoupled, as required. For the Gaussian model, 
there is no separate spine component, so the inner region 
has a constant velocity (5 = j3i everywhere (i.e. f3 s = 1). An 
unphysical acceleration is required in the shear layer at the 
flaring point in both classes of model: this is an inevitable 
consequence of the increase in sidedness ratio. We discuss 



this problem and a possible solution in Section 5.1 



parameters required to specify the functional form of the 
variation. Table ^ summarizes the details, as follows: 

Column 1: Symbol (as defined in the text), 

Columns 2-5: the values of the quantity at the four stan- 
dard locations (blank if not used), 

Column 6: Any other parameters needed to quantify the 
variation. 

Columns 7-9: The functional forms of the variation in 
the three regions. 

Similarly, the functional forms used to describe transverse 
variations in the shear layer are listed in Table |B| There are 
no variations across the spine in SSL models. 



3.4 Velocity field 

We have chosen to model the velocity field as a separable 
function (3(p,s) = f3 p (p)/3 s (s) with /3 S (0) = 1. The inner 
region is faint and poorly resolved so, in the absence of ev- 
idence to the contrary, we assume that the on-axis velocity 
is constant there. To generate the sideness profile of the rest 
of the jet (Fig. ^|), the on-axis velocity must remain fairly 
constant throughout much of the flaring region, drop rapidly 
just before the outer boundary and then fall smoothly and 



3.5 Emissivity 

As with the velocity, we use a separable function for the 
rest-frame emissivity: e(p, s) = e p (p)e s (s). We found that 
very different gradients of the on-axis emissivity e p (p) were 
required in the three regions, and therefore adopted a power- 
law form, with different exponents allowed for the regions, 
and for the spine and shear layer (Table |^). One additional 
parameter is needed to set the relative emissivity of spine 
and shear layer at a fiducial point. We enforce continuity 
at the outer boundary, but could not fit the data without 
intro ducin g a discontinuity at the inner boundary (see Sec- 
tion 1.3.5 ). The transverse variation has the same form as 
that assumed for the velocity: constant in the spine, with 
a linear (SSL) or truncated Gaussian decrease in the shear 
layer to a fraction e(p) at the jet edge. The absolute value of 
the emissivity is determined by normalizing to the observed 
flux density. 



3.6 Field ordering 

Clues to the three-dimensional structure of the magnetic 
field come from the differences in polarization betw een the 
main and counter-jets, as summarize d in Section 2.4|. We ini- 



tially tried the structure proposed by Laing (1995 ) in which 
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Table 4. Summary of the functional variations of velocity, emissivity and field ordering parameters along the model jets. 



Quantity 



o-n 



Free parameters 
r\ ro rf Other 



Inner 



Functional dependences 
Flaring Outer 



Hp) 
Hp) 



hip) 

k p {p) 



Velocities 
ft 0i 



ft) ft 

t'n 



H 



ft 



Emissivity (for each of spine and shear layer) 

g, E h Et, E g(p/n)- E ' 
ei e 1 



bo + b ± p 
vi 



■ b 2 p l 



ei 



Field component ratios (for each of spine and shear layer) 
3i h 30 if ji h - 

k\ ki ko kf ki ki 



(r— ri)C«a- 




r 0~ r l 








(l — ri)(e - 


-ei) 


r -ri 




(r-ri)(jo- 


Ji) 


ro-rx 




(l — ri)(fe - 


-fei) 



c exp(-cip) 

(r /r 1 )^ E !(p/r )" 
eo 



Jo - 



(r-i-o)(j f -jo) 

r-f-ro 
(r-r )(fc f -fc ) 



a fast, transverse-field spine with equal radial and toroidal 
field components is surrounded by a slower longitudinal-field 
shear layer. This might naively be expected from a combina- 
tion of expansion and interaction with the external medium. 
Such a structure produces a transverse apparent field on-axis 
and a longitudinal field at the edge, together with a transi- 
tion from longitudinal to transverse apparent field on-axis, 
both as observed (Fig. ^|). It could be rejected for 3C31, 
however, because it always requires the transition from lon- 
gitudinal to transverse apparent field to occur closer to the 
nucleus in the main jet (where Doppler boosting makes the 
faster spine emission with its transverse apparent field rela- 
tively more prominent). The opposite is observed. The high 
degree of polarization observed in the outer counter-jet is 
also inconsistent with such a field configuration. 

This suggested a model in which both the spine and the 
shear layer have toroidal and longitudinal field components 
(of roughly equal magnitude) but t he radial co mponent is 
everywhere very small (model B of Laing 1981). The field 



Table 5. Summary of the functional variations of velocity, emis- 
sivity and field ordering parameters across the model shear layers. 



is then two-dimensional, in sheets wrapped around the jet 
axis. The apparent field is always longitudinal (with the the- 
oretical maximum degree of polarization, po) at the edges 
of the jets, but can be either longitudinal or transverse on 
the axis, depending on the relative magnitudes of the two 
components and the angle to the line of sight. If this an- 
gle and the flow velocities are adjusted appropriately, then 
aberration can act so that the field sheets are seen face-on 
in their rest frames in the main jet (giving a low degree 
of polarization), but side-on in the counter-jet (leading to 
high polarization and a transverse apparent field). Models 
of this type produce much more realistic polarization distri- 
butions, especially when the ratio of toroidal to longitudinal 
field increases with distance from the nucleus, but still fail 
in two important respects. First, the field transition region 
in the main jet is too close to the nucleus and too short. Sec- 
ond, a high degree of polarization is predicted at the edge 
of the flarin g re gion, where the observed values are quite 
low (Section 2.4). The solution to both problems is to allow 
a radial field component which increases from zero close to 
the axis to a finite value at the edge of the jet. This edge 



Quantity 



Model Functional variation 



Velocity 

s (s)(p>pi) SSL 1 + 

Gaussian exp[— s 2 ln«(p)] 
ft(s) (p < pi) SSL Vi 

Gaussian 1 



.(«) 



3s(s) 



Emissivity 
SSL 1 + [e(p) - l]a 

Gaussian exp[— s' 2 In e(p)] 

Radial/toroidal field ratio 



value must vary along the jet in such a way that the field is 
essentially isotropic at the boundary in parts of the flaring 
region, but the radial component vanishes at large distances 
from the nucleus. In contrast, we found no evidence for any 
transverse variation of the longitudinal/toroidal ratio in the 
shear layer. 

The functional forms are again given in Tables ^ and 
|E| We use the ratios of rms field components j(p, s) (ra- 
dial/toroidal) and k(p,s) — k p (p) (longitudinal/toroidal), 
with no transverse variation in the spine for SSL models. 
We chose j(p,s) = j P (p)js(s) with j s (s) = s p for the ra- 
dial/toroidal ratio. If the functional forms given for the outer 
region are negative, the corresponding values of j or k are 
set to zero. 



3.7 Model integration 

The principal steps in calculating the brightness distribu- 
tions are as follows: 
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(i) Construct grids to match the observations at each of 
the two resolutions. 

(ii) At each grid point, determine whether the line of sight 
passes through the jet. If so, calculate the integration limits 
corresponding to the outer surface of the jet and, if relevant, 
the spine/shear-layer interface. Separate ranges of integra- 
tion are required to avoid discontinuities in the integrand. 

(iii) Integrate to get the Stokes parameters /, Q and U 
using Romberg integration. The steps needed to determine 
the integrand are outlined below. 

(iv) Add in the core as a point source. 

(v) Convolve with a Gaussian beam to match the resolu- 
tion of the observations. 

(vi) Evaluate x 2 over defined areas, using an estimate of 
the "noise level" derived as described later. 

In order to determine the I, Q and U emissivities at a point 
on the line we follow an appro ach described in detail in [Laing 



(2002) and based on that of Matthews fc Scheuer (1990) 



We neglect synchrotron losses, on the grounds that the ob- 
served spectrum is a power law with a = 0.55 between 1.4 



and 8.4 GHz (and e xtends to much higher frequencies; Hard- 
castle et al. 2002). The emissivity function e oc noB L+a 



where B is the total field and no is the normalizing cons tant 



in the electron energy distribution as defined in Section 3.1 



The observed emissivity can be calculated in the formalism 
developed by Laing (2002) by considering an element of 



fluid which was initially a cube containing isotropic field, 
but which has been deformed into a cuboid by stretching 
along the three coordinate directions by amounts propor- 
tional to the field component ratios in such a way that the 
value of e is preserved. We calculate the synchrotron emis- 
sion along the line of sight in the fluid rest frame, thus taking 
account of aberration. 

The main steps in the calculation are: 

(i) Determine coordinates in a frame fixed in the jet, in 
particular the radial coordinate p and the streamline index 
s, numerically if necessary. 

(ii) Evaluate the velocity at that point, together with the 
components of unit vectors along the streamline coordinate 
directions (and hence the angle between the flow direction 
and the line of sight ip). Derive the Doppler factor D — 
[r(l — /Jcosi/')]" 1 and hence the rotation due to aberration 
(sin?// = Dsmtp, where ip' is measured in the rest frame 
of the jet material). Rotate the unit vectors by ip — ip' and 
compute their direction cosines in observed coordinates. 

(iii) Evaluate the emissivity function e and the rms com- 
ponents of the magnetic field along the streamline coordi- 
nate directions (normalized by the total field). Scale the 
direction cosines derived in the previous step by the cor- 
responding field components, which are + j 2 + fc 2 ) 1 ' 2 
(radial), 1/(1 + j 2 + k 2 ) 1/2 (toroidal) and k/(l + j 2 + k 2 ) 1/2 
(longitudinal) in the notation of the previous section. 

(iv) Evaluate the position angle of polarization, and the 
rms field components along the major and minor axes of 
the probability d ensity functio n of the field projected on the 
plane of the sky flLaing 20021) . Multiply by e(p, s)D 2+a , to 
scale the emissivity and account for Doppler beaming. 

(v) Derive the tot al and polariz ed emissivities using the 
expressions given by Laing (2002) and convert to observed 
Stokes Q and U. 



3.8 Fitting and optimization 

Our basic approach is to minimize x 2 between the model 
predictions and the data, summing values for the three in- 
dependent Stokes parameters, /, Q and U. The value of the 
"noise" on the observed images is important in the opti- 
mization process, as sums of x 2 over different areas need 
to be added with the appropriate weights to ensure that 
the data are fitted sensibly over the full range of resolutions 
available. The "noise" is dominated by small-scale intensity 
fluctuations - knots and filaments - whose amplitude is un- 
known a priori. Our best guess at their level comes from a 
measure of the deviation from axisymmetry. The "noise", 
E, is taken to be l/\/2 times the rms difference between the 
image and a copy of itself reflected across the jet axis. This 
is always much larger than the off-source rms. Any contri- 
bution from deconvolution artefacts will also be included in 
this estimate. Some components of the small-scale structure 
will result in mirror-symmetric features in the brightness 
distribution (e.g. the bright arc in the main jet; Fig. hi), and 
we will therefore underestimate E. 

We fit to the 0.25-arcsec images over the area covered 
by the model from 0.5 to 4.1 arcsec from the core. This ex- 
cludes the core itself, and covers all of the area where signif- 
icant polarized emission is detected at this resolution. Fur- 
ther out, the signal-to-noise ratio for these images (especially 
in linear polarization) is too low to provide an effective con- 
straint, so we fit to the 0.75-arcsec images. Fits made using 
the high-resolution images alone are consistent with those 
that we describe here, but are less well constrained. \ 2 IS 
computed only over the model area and we evaluate it at at 
grid-points chosen to ensure that all values are independent. 
There are 1346 independent points, each with measurements 
in 3 Stokes parameters. Of these, 44, 162 and 1140 are in 
the inner, flaring and outer regions, respectively. 

We have optimized the models over the whole area and 
with one or more of the brightest small-scale features ex- 
cluded from the \ 2 calculation. The derived parameters did 
not vary by appreciable amounts, but exclusion of the ob- 
vious "arc" in the main jet (Fig. [j]) somewhat reduced the 
final x 2 - Given that we are effectively averaging over many 
small-scale filaments in the brightness distribution, we have 
no physical reason to remove the brightest few, but it is 
reassuring that the results are insensitive to their exclusion. 

In order to optimize the model parameters, we use the 



down hill simplex method of Nelder & Mead (Press et al 



1992). This usually converges in 150 - 200 iterations given 



reasonable starting parameters. 



3.9 Uniqueness 

As with any model-fitting procedure, questions of unique- 
ness must be considered. Our approach is an advance on 
previous attempts at jet velocity estimation in several re- 
spects: 

(i) We seek to fit a large quantity of well-resolved two- 
dimensional data, rather than one-dimensional profiles. 

(ii) We have detected both jets at all distances from the 
nucleus, and do not have to cope with upper limits. 

(iii) We fit linear polarization (Stokes Q and U) and to- 
tal intensity (Stokes I) simultaneously with a small number 
of free model parameters. Although this introduces further 
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Figure 5. Grey-scale images of \ 2 = (Smodel ~ •S'obsorvod) 2 /^ 2 , where E is the "noise level" defined in the text. Left: images with a 
resolution of 0.75arcsec FWHM covering an area of ±27arcsec from the nucleus. Right: images with 0.25arcscc FWHM of the inner 
±7.5arcsec. The labelled bars indicate the grey-scale ranges: these are different for the two resolutions. The boxes show the areas over 
which x 2 was summed at each resolution to assess goodness of fit. From the top: x 2 images for Stokes /, Q and U; total intensity for 
the same area. \ 2 values are not plotted for Stokes I in the immediate vicinity of the core. 



degrees of freedom in order to describe the field anisotropy, 
we find that the form of the jet velocity field is as severely 
constrained by the observed polarization data as by the jet 
sidedness - the more traditional quantity used to infer jet 
velocities. 



Model images that even qualitatively resemble the observa- 
tions are hard to find. Although the downhill simplex algo- 
rithm is not guaranteed to converge on a global minimum in 
X 2 , we experimented with a wide range of initial conditions 
and found no other significant minima. We are therefore con- 
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Gaussian model 




SSL model 




Figure 6. Contours of total intensity at a resolution of 
0.75arcsec, covering ±27arcsec from the nucleus. The contour 
levels are: -1, 1, 2, 4, 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96, 
104 X 40 /ijy/beam area. From top: model with Gaussian profile, 
model with spine/shear layer, VLA data. 



Table 6. Summary of reduced x 2 values. 



Model type 


No blanking 


Arc blanked 


SSL 


1.71 


1.51 


Gaussian 


1.80 


1.60 



fident that the parameters given in the next section describe 
a unique solution. 

4 MODEL RESULTS 

4.1 Comparison between models and data 

4.1.1 \ 2 values 



Gaussian model 




SSL model 




1 arcsec 



Figure 7. Contours of total intensity at a resolution of 
0.25 arcsec. The plot covers ±10 arcsec from the nucleus. The 
contour levels are -1, 1, 2, 3, 4, 6, 8, 10, 12, 14, 16, 20, 24 X 
30 /^Jy/beam area. From top: model with Gaussian profile, model 
with spine/shear layer, VLA data. 

We have optimized models with spine/shear layer and 
Gaussian transverse profiles. The resulting reduced \ 2 val- 
ues, with and without blanking of the brightest arc in the 
main jet, are given in Table |^ and images of \ 2 axe shown 
in Fig. |^. Both models fit the large-scale total intensity and 
polarization distributions well. Given that the "noise levels" 
are crudely estimated (and likely to be too low), and that 
the "noise" shows large-scale correlation with clearly non- 
Gaussian statistics, it is not unexpected that the reduced 
X 2 « 1.5-1.8 is inconsistent with a formal fit. 

In what follows, we concentrate on the SSL model as the 
best description of velocity, emissivity and field ordering re- 
gardless of the underlying physics. The Gaussian equivalent 
has a smaller number of free parameters (Table [?]). It gives 
a slightly, but significantly worse fit, except in the inner re- 
gion, where it fails seriously (albeit with little effect on the 
overall \ 2 )- 
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Figure 8. Total intensity profile along the jet ridge line at 
0.75 arcsec resolution. Full line: data; dashed line: SSL model; 
dotted line: Gaussian model. 



CO 
CD 

CO 

6 

CO 

<u 



0.8 



0.6 



0.4 




0.2 





-4 -2 2 4 
Distance from nucleus (arcsec) 

Figure 9. Total intensity profile along the jet ridge line at 
0.25 arcsec resolution. Full line: data; dashed line: SSL model; 
dotted line: Gaussian model. 




Figure 10. Grey-scale images of the observed (bottom) and SSL 
model (top) sidedness ratios at a resolution of 0.75 arcsec. Both 
images use the same transfer function, which has been optimized 
to emphasize the variations in sidedness in the outer part of the 
jet. The sidedness images for the SSL and Gaussian models are 
very similar, so the latter is not shown. 



4-1.2 Total intensity images and profiles 

Figs | and @ show the predicted and observed images of 
Stokes / at 0.75 and 0.25 arcsec resolution, respectively. Pro- 
files of I along the jet axis are given at these resolutions in 
Figs |i| and []. The differences between the main and counter- 
jets are emphasized in sidedness images and profiles (Figs |lo| 
and0). 



4-1.3 Fitted total intensity features 

The following features of 3C31 can be accurately reproduced 
by our chosen fitting functions after optimization: 

(i) Both jets are initially faint and brighten at the be- 
ginning of the flaring region, where significant deceleration 
begins. 

(ii) The brighter jet has a more centrally-peaked bright- 
ness distribution, while that of the counter-jet is much flat- 
ter. 

(iii) The jets become more equal in brig htness further 
from the nucleus as they decelerate (Fig. |l0|) . 

(iv) The on-axis sidedness ratio remains high («13) over 
most of the flaring region, and drops abruptly at 5 arcsec 
from the nucleus. Thereafter, it declines slowly and mono- 
tonically but the main jet remains appreciably brighter than 
the counter-jet on-axis (Fig. pd| ). 



3C31 jet model 13 



25 



CO 

u 



* 20 
+-) 

W 
S=! 
CD 

fl 15 



CD 

T 

Sh 

CD 

p 

o 
o 



CD 

!-3 



10 







I I I I I I I I I I I I I I I I I I I I I I I I I I 




5 10 15 20 25 
Distance from nucleus (arcsec) 

Figure 11. Profiles of jet/counter-jet sidedness ratio along 
the axis at a resolution of 0.75 arcsec. Full line: data; dashed: 
spine/shear layer model; dotted: Gaussian model. 



The differences between the spine/shear layer and 
Gaussian models are at a low level. The former allows a 
lower emissivity in the spine, which leads to a flatter trans- 
verse intensity profile that agrees better with the data. 



4-1-4 Polarization images and profiles 

Fig. ^ shows the predicted and observed degrees of polar- 
ization p — (Q 2 + U 2 ) 1 / 2 /I at both resolutions, with the I 
images below them for reference. The degree and direction 
of polarization are represented in Figs |l|and[lj by vectors 
whose magnitudes are proportional to the degree of polar- 
ization and whose directions are those of the apparent mag- 
netic field (i.e. rotated by 90° from the E-vector direction, 
after correcting the observations for Faraday rotation). Lon- 
gitudinal and representative transverse profiles at the lower 
resolution are displayed in Figs |l5| and respectively. 



4-1.5 Fitted polarization features 

Our choice of fitting functions also reproduces the following 
features of the polarization distribution: 

(i) There is a V-shaped region at the onset of the flaring 
of the main jet in which the polarization is close to zero on- 
axis and rises to values «40% with longitudinal apparent 
field at the edge (Fig. 

(ii) Between 4 and 8 arcsec from the nucleus in the flaring 
region, the main jet polarization drops to a very low level 
over the entire width (Fig. 

(iii) In contrast, the corresponding region of the counter- 
jet shows a transverse apparent field with monotonically in- 



creasing polarization on-axis with a much lower degree of 
polarization towards the edge (Fig. |l6| ). 

(iv) In the outer region, the degree of polarization in- 
creases monotonically with distance from the nucleus in both 
jets, with a transverse apparent field, but the degree of po- 
larization is always higher in the counter-jet (Fig. p"5| ). 

(v) Both jets show longitudinal apparent field at their 
edges, with a degree of polarization approaching 70% at the 
extreme outer edge (Figs ^ and ^J) . 

4-1.6 Features that cannot be fitted well 

The models are in principle incapable of fitting non- 
axisymmetric or small-scale features. The most important 
examples of these, emphasized in the plots of \ 2 f° r Stokes 
/ and U (Fig. ^) are as follows: 

(i) The inner and flaring regions of the main and counter- 
jets have fine structure consisting of small numbers of dis- 
crete knots. These are modelled as continuous features with 
the correct mean level (Fig. ^). 

(ii) The observed apparent magnetic field direction is 
oblique to the jet axis in the centre of the flaring region 
of the main jet (Fig. |l4| ): this cannot be reproduced in any 
purely axisymmetric model. 

(iii) The prominent arc of emission 20 to 24 arcsec from 
the nucleus in the main jet is not reproduced either in total 
intensity or linear polarization. 

In addition, there are small but significant deviations 
between observed and modelled polarization patterns in the 
flaring region: 

(i) The apparent field vectors between 5 and 7 arcsec from 
the nucleus in the main jet diverge more from the axis than 
is predicted (Fig. |l4|). The degree of polarization also ap- 
pears to be underestimated, but the signal-to-noise ratio at 
0.25 arcsec FWHM is quite low and that the observed vectors 
are blanked on polarized flux. There is therefore a tendency 
for the degree of polarization to be spuriously high for the 
plotted vect ors. This effect has been corrected to first order 
(Section 2.3), but some residual remains. 

(ii) The degree of polarization along the ridge-line of the 
countcr-jet is underestimated, sig nificantly so between 3 and 
10 arcsec from the core (Figs [L5| and |l6| ). 



4.2 Model parameters and confidence limits 

The optimization problem is complicated, and estimates of 
some of the parameters are strongly correlated. In addition, 
we do not know the statistics (or even the rms level) of the 
"noise" a priori and we have imposed additional constraints 
by our choice of fitting functions. The \ 2 statistic is effective 
in optimizing the fit, but assessing confidence limits (e.g. by 
a Bayesian analysis or using bootstrap techniques) would be 
far from straightforward. We have instead adopted a simple 
ad hoc procedure, by which we vary a single parameter until 
the fractional increase in the \ 2 values for I or Q and U in 
one of the inner, flaring or outer regions corresponds to the 
formal 99% confidence limit for independent Gaussian errors 
and the appropriate number of degrees of freedom. Most pa- 
rameters affect the fit significantly only for part of the jet, 
or for a subset of the Stokes parameters, so this approach 
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Figure 12. Grey-scale images of the total intensity and degree of polarization, p, at a resolution of 0.75arcsec (left) and 0.25arcsec 
(right). The labelled bars give the grey-scale levels for p and are different for the two resolutions. From top: polarization data, polarization 
of model with spine and shear layer, total intensity data. The data and models have been blanked wherever the polarized signal is < 3<r 
or the total intensity is < 5<r, using the values of off-source noise, <r, given in Table bl 



is superior to one based on the total \ 2 ■ The estimates are 
qualitatively reasonable, in the sense that varying a param- 
eter by its assigned error leads to a visibly unacceptable fit, 
and we believe that they give a good general impression of 
the range of allowed models. The numerical confidence levels 
should not be taken too seriously, however. 

Table gives the fitted parameters and error estimates 
for the spine/shear layer and Gaussian models. The param- 
eters are the angle to the line of sight, the spine opening 
angles and those defined in Tables W - 0. The columns are: 



(iii) The minimum value allowed by the \ 2 recipe given 
earlier, again for the SSL model. 

(iv) The corresponding maximum value. 

(v) The best fit for the Gaussian model (the allowed 
ranges are very similar to those for the SSL model). 

In general, the parameters for the Gaussian and SSL models 
are very similar and always agree within the quoted errors 
- the contribution of the spine component to the emission 
(and therefore to the \ 2 value) is quite small. 



(i) A description of the parameter. The symbols are those 
used in Tables § - |[ 

(ii) The best fit for the SSL model 
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Table 7. Fitted parameters and error estimates. 



Quantity SSL Gauss 

opt minP maxr 



Angle to line of sight 8 (degrees) 


52.4 


48.9 


54.1 


51.4 


Jet half-opening angles (degrees )p| 










inner region 


6.7 






6.6 


outer region £ 


13.2 






13.0 


Boundary positions (kpc) 










inner r\ 


1.1 






1.1 


outer T'o 


3.5 






3.6 


arbitrary fiducial Tr 


9.6 






9.8 


Spine half opening angles (degrees) 










inner region £j 


4.06 


3.1 


6 5 




outer region C^o 


2.79 


0.7 


1.5 




On — axis velocities / c 










inner jet /3; 


0.87 


0.83 


0.93 


0.20 


inner boundary [3\ 


0.77 


0.68 


0.83 


0.76 


outer boundary [3q 


0.55 


0.45 


0.63 


0.54 


outer fiducial /3f 


0.28 


0.25 


0.33 


0.27 


velocity exponent H 


9.5 


3.6 


> 


8.8 


Fractional velocity at edge of jetQ 










inner jet u; 


0.06 


0.0 


1.15 




inner boundary v\ 


n 7 a 


n i 


l.OV 


n Q7 


outer boundary vq 


U.O f 


U.ol 


U.o / 


U.Do 


On — axis emissivity exponents 










inner spine _E; 


1.96 


< 


2.3 




flaring spine E{ 


2.52 


1.9 


2.9 




outer spine E a 


2.14 


1.1 


3.8 




inner shear layer E\ 


1.33 


< 


2.2 


0.75 


flaring shear layer Ef 


3.10 


2.9 


3.4 


3.08 


outer shear layer E Q 


1.42 


1.33 


1.54 


1.44 


Fractional emissivity at edge of jet 










inner boundary e\ 


0.27 


0.05 


0.52 


0.37 


outer boundary eo 


0.20 


0.09 


0.28 


0.26 


Shear layer / spine emissivity 


2.11 


1.5 


3.1 




ratio at inner boundary 










Emissivity ratio at inner boundary 










(inner / flaring region) 










spine g 


0.37 


0.13 


0.53 




shear layer g 


0.04 


0.003 


0.08 


0.05 



a The symbol < means that any value smaller than the quoted maximum is allowed. 

^ The symbol > means that any value larger than the quoted minimum is allowed. 

c Opening angles and boundary locations are given in the jet coordinate system. The jet 

opening angles and the boundary locations are determined by the outer isophotes once the 

angle to the line of sight is specified , so no errors are quoted. 

^ The upper limits on the fractional velocity at the edge of the jet in the inner region and 
at the inner boundary are set not by the \ 2 constraint but rather by the condition that the 
velocity must be < c. 
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Figure 13. Vectors whose lengths are proportional to the degree of polarization, p, and whose directions are those of the apparent 
magnetic field superimposed on grey-scale images of total intensity at 0.75 arcsec resolution. Top: spine/shear-layer model; bottom: VLA 
data. The vectors are blanked as in Fig. 



4.3 Model description 

In this sub-section, all of the plots represent a plane con- 
taining the jet axis (not projected on the sky) and distances 
are given in linear units. 



4-3.1 Angle to the line of sight 

The angle to the line of sight is one of the few parameters 
that affects the total and linearly-polarized emission on all 
scales, and is therefore particularly well constrained. The 
best value for both Gaussian and SSL models is 9 w 52°, 
with uncertainties of a few degrees. 



4.3.2 Velocity field 

The inferred velocity field is shown as contour plots in 
Fig. [Tt] and as longitudinal profiles in Fig. The on-axis 
velocity of the inner jet is poorly constrained, although ex- 
tremely low values (/3 < 0.4) are ruled out. The best fits 
show an abrupt decreas e in velocity across the inner bound- 
ary (Fig [l^; see Section 



5.1 



for a more detailed discussion), 
but continuity cannot be entirely ruled out. Further out, 
the velocity field is much better constrained (±0.1c and 
±0.05c, respectively, for the flaring and outer regions) and 



both models agree almost exactly. The on-axis velocity re- 
mains roughly constant {13 ~ 0.77) between 1 and 2.5 kpc 
and then drops abruptly to f3 ss 0.55 at the outer boundary. 
Quantitatively, the exponent H in the velocity law for the 
flaring region (Table W) is required to be > 3. Thereafter, j3 
declines smoothly to ~0.22 at 12 kpc. 

The transverse velocity profile is hardly constrained at 
all in the inner region. The best fits in the flaring and outer 
regions require an edge velocity close to 0.7 of the central 
value (Figs [l7| and Jig), independent of distance from the 
nucleus. The error analysis shows, however, that we cannot 
exclude a flat-topped profile at the inner boundary, so some 
evolution of the profile along the jet could occur. Very low 
velocities at the edge of the jet are not consistent with the 
observed sidedness ratios in these regions. 



4-3.3 Emissivity 

The spatial variation of noB 1+a (proportional to the rest- 
frame emissivity) is shown as a grey-scale image in Fig. |l^ 
and as longitudinal profiles in Figs The emissivity is 
again poorly constrained in the inner region, and we can 
only exclude very steep decreases with distance. For the 
shear layer, the dependence of emissivity on distance is very 
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Table 7. Fitted parameters and error estimates (continued). 



Quantity SSL Gauss 

opt min max 



RMS field ratios (shear layer) 



radial/toroidal 



inner jet centre j\ 


0.37 


0.0 


> 


0.38 




0.0 






0.0 


inner boundary centre j± 


0.93 


0.3 


1.4 


0.78 


inner boundary edge j 


0.0 






0.00 


outer boundary centre jo 


1.00 


0.52 


1.38 


0.92 


outer boundary edge 


0.0 






0.00 


fiducial distance centre jf 


0.0 


0.0 


0.62 


0.24 


fiducial distance edge 


0.0 






0.00 


index p 


0.53 


0.3 


1.5 


0.41 


longitudinal / toroidal 










inner jet fc; 


1.23 


0.2 


2.3 


1.43 


inner boundary k\ 


1.16 


1.05 


1.35 


1.17 


outer boundary ko 


0.73 


0.63 


0.80 


0.82 


fiducial distance kf 


0.50 


0.41 


0.58 


0.54 


RMS field ratios (spine) 










radial / toroida^] 










inner jet j; 


0.0 


0.0 


1.5 




inner boundary jr'i 


0.0 


0.0 


1.3 




outer boundary jo 


0.0 


0.0 


1.9 




fiducial distance j{ 


0.0 


0.0 


10.0 




longitudinal / toroidal 










inner jet fc; 


1.75 


1.1 


2.1 




inner boundary ki 


1.06 


0.7 


1.8 




outer boundary ko 


1.40 


0.8 


4.0 




fiducial distance fcf 


0.84 


0.0 


8.0 





e The radial/toroidal ratios always vary from at the spine/shear-layer interface (SSL) or 
axis (Gaussian) to a maximum value at the edge of the jet (Table ^).The values quoted are 
for the edge and centre of the shear layer. 

f radial / toroidal field ratios for the spine in the SSL models are consistent with but poorly 
constrained everywhere (to the extent of having negligible influence on the x 2 values). The 
relevant parameters were fixed at throughout the optimization process. 



different for the flaring (oc p -3 ' 1 ) and outer (oc p~ 1A ) re- 
gions. The relative contribution from the spine is small, so 
its emissivity could vary either as a single power-law or in 
a manner closer to that of the shear layer. The fractional 
emissivity at the edge of the shear layer is consistent with a 
constant value of «0.2 in the flaring and outer regions, al- 
though a much wider range (including zero) is allowed at the 
inner boundary. There is strong evidence for a discontinuity 
in emissivity at the inner boundary for the shear layer. 



4-3-4 Magnetic-field structure 

Very little polarized flux is either predicted or observed to 
come from the inner region, and no polarized emission is 
detected from the inner counter-jet. For this reason, the field 
ratios are essentially unconstrained there (Table Q). 

The toroidal field component is the largest over most of 



the flaring and outer regions, increasing from «0.6 to ~0.9 
of the total between the inner boundary and the end of the 
modelled region (Fig. |2l| ). The longitudinal component, con- 
versely, decreases from ~0.7 to «0.4 over the same distance. 
We found no evidence for any variation of the longitudi- 
nal/toroidal field ratio across the shear layer. The compo- 
nent ratios for the spine are not well determined (Table 
and could quite plausibly be identical to those for the shear 
layer. 

Close to the edge of the jet in the flaring region, the ra- 
dial component becomes appreciable, reaching a maximum 
of «0.6 at the edge on the outer boundary (marked by the 
arrows in Fig. |2lj). Thereafter, it decreases rapidly with dis- 
tance along the jet, becoming negligible by 9.5 kpc from the 
nucleus in the SSL model. The radial component, unlike the 
other two, increases with radius (Fig. Ell). As mentioned 
earlier, this variation is required in order to achieve even 
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SSL model 




0.6 



Data 



p = 1 




5 arcsec 



Figure 14. Vectors whose lengths are proportional to the degree 
of polarization, p, and whose directions are those of the appar- 
ent magnetic field superimposed on contours of total intensity 
at 0.25 arcsec resolution. Top: spine/shear-layer model; bottom: 
VLA data. The observed vectors are blanked on polarized and 
total intensity, as in Fig. |l2| but the model vectors are plotted 
wherever the total intensity exceeds 5cr, regardless of polarized 
flux. 



a qualitative fit to the observed polarization in the flaring 
region. 



5 DISCUSSION 

5.1 The inner region 

Our estimates of the angle to the line of sight and the jet 
ve locity in the inner r egion are entirely consistent with those 
of Lara et al. (1997) for the jets on parsec scales. For our 



preferred value of 6, their velocity range is 0.81 < (3 < 0.998, 
consistent with our central velocity for the inner region (/3 ss 
0.87), but also allowing significant deceleration from parsec 
to kiloparsec scales. 

The inner region poses a problem for any decelerating- 
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Figure 15. Profile of the degree of polarization along the jet 
ridge line at 0.75 arcsec resolution. Full line: data; dashed line: 
spine/shear layer fit; dotted line: Gaussian fit. 

jet model, in which the jet-to-counterjet intensity ratio (sid- 
edness) must decrease with distance from the nucleus. We 
would expect the sidedness ratio to have a maximum in the 
inner region, in which case the counter-jet would be invisi- 
ble. In fact, the inner region clearly has a lower sidedness, 
on average, than that in the flaring region. The brightness 
distribution is dominated by a few knots, so one possibility 
is that we are being fooled by stochastic variations. Alter- 
natively, there could be a small amount of relatively slow- 
moving material in the shear layer, surrounding a fast spine. 

If a very slow component also exists further out, it can- 
not have a noticeable effect on the brightness distribution: 
the fact that the sidedness ratio at the edges of the jet in 
the flaring region differs significantly from unity requires 
that the emissivity of any very slow component becomes 
insignificant on large scales. This component therefore has 
a negligible effect on the fits beyond the inner jet, and its 
properties are constrained only by the intensity fits in the 
inner region. The faster component of the flow that dom- 
inates the outer jets is relatively faint in the inner region 
(because its Doppler factor and intrinsic emissivity are both 
lower than on larger scales), so the modelling of the inner 
region is almost decoupled from that of the rest of the jets 
in the spine/shear- layer fits. 

In the best fit SSL model, the slow component has been 
introduced by assigning a velocity of 0.06c to the shear layer 
while the spine has a velocity of 0.87c, i.e. we have allowed 
the slow component in the inner region to substitute for the 
faster-moving shear layer component that is required to ex- 
plain the flaring region via an unphysical jump condition at 
the boundary. We emphasize that higher resolution imag- 
ing of the inner jets is needed to obtain firmer constraints 
on the transverse velocity distribution in this region, and to 
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Figure 16. Example transverse profiles of the degree of polar- 
ization at distances of 8 arcsec from the core in the flaring region 
(left panels) and f 6 arcsec from the core in the outer region (right 
panels). Top panels: main jet; bottom panel: counter-jet. Full 
line: data; dashed line: spine/shear-layer fit: dotted line: Gaus- 
sian fit. The profiles are blanked as in Fig. n2] and the resolution 
is 0.75 arcsec. 



explore how this distribution evolves as the jet enters the 
flaring region. Our best guess at the geometry is sketched 
in Fig. hi If this picture is correct, there must still be an 
increase in emissivity at the flaring point, but the values of 
g quoted in Table Q will be inaccurate. Finally, we note that 
an additional component with f3 significantly higher than 0.9 
would be severely Doppler-dimmed even in the approaching 
jet, and therefore very difficult to detect. 



5.2 The onset of flaring and deceleration 

We have shown that the onset of deceleration is marked by 
a large increase in rest-frame emissivity and a major change 
in the jet collimation. It is not merely that the jet becomes 
gradually brighter as it decelerates and Doppler suppression 
is reduced: there is also a discontinuity at the inner bound- 
ary. 

One possibility is that the jet is supersonic, over- 
pressured and expanding freely in the inner region. In that 
case, the internal pressure would fall until it drops below 
that of the extern al medium, at which point a reconfine- 



ment shock forms (Sanders 1983). The reconfinement shock 
is followed by a second shock at which the jet becomes over- 
pressured with respect to the external medium and it this 
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Figure 17. Contours of the velocity field for the fitted models. 
Top: model with Gaussian profile; bottom: model with spine and 
shear layer. Contours are shown at intervals of 0.05c and fiducial 
contours are labelled in the outer jets. The panels correspond to 
the same area projected on the plane of the sky and differ very 
slightly in size because the values of 9 are not identical in the two 
models. Both cover very nearly f 2 kpc in the plane of the jet. 
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Figure 18. Profiles of the velocity along streamlines for the SSL 
model. Full line: spine (on-axis) and shear layer s = 0; short dash: 
shear layer s = 0.5; long dash: shear layer s = f (jet edge). 
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Figure 19. Grey-scale image of \og(noB 1+a ), derived from the 
model emissivity, for no in m -3 and B in T. Top: model with 
Gaussian profile; bottom: model with spine/shear layer. The areas 
covered by the two plots are not quite the same because the two 
models have different angles to the line of sight, but the linear 
scales are identical. 



feature which is most plau sibly identified with the flaring 
point. For a relativistic jet, Komissarov (1994) shows that 
the shock forms at a distance 



•2-shock 



/ 2$ 



\3TTPcxtC 



1/2 



where $ is the energy flux throu gh the jet and p cxt ~ 
3x 1CP 11 Pa (Hardcastle et al. 2002) is the external pressure. 
This would be consistent with the observed inner boundary 
distance of 1.1 kpc for an energy flux of « 5 x 10 37 W, some- 
wh at higher than that th e value of ~ 1 x 10 37 W estimated 
by Laing & Bridle (2002) from a conservation- law analysis. 

We see no evidence for any simple shock structure at 
the inner boundary, although the emission there is not com- 
pletely resolved and there are obvious (non-axisymmetric) 
knots at the beginning of the flaring region. If the inner re- 
gion is in free expansion, we can estimate the initial Mach 
number of the flow from the opening angle: arctan(£i) « A4 
where M = (T0) / (T 3 3 ) is the generalized Mach number 
defined by Konigl (198C ), 0., = c.,/c, c. s is the internal sound 
speed, < ,nd F 3 = ( 1 — 0^ l^ 1 ^ 2 . The observed value of ^ = 6.7° 



corresp mds to M ~ 8.5 and hence to T ~ 6.1 if the inner jet 



has the sound speed c s = c/vo of an ultra-relativistic gas. 
This initial velocity is considerably faster than we have in- 
ferred for the inner region but, as mentioned in Section |5.l[ 
we cannot exclude the presence of such higher- velocity ma- 
terial there. 
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Figure 20. Profiles of tlqB +a inferred from the rest-frame emis- 
sivity for two streamlines of the SSL model compared with those 
predicted by adiabatic models (no and B are in m — 3 and T, re- 
spectively). Full lines: SSL; long dashes: parallel-field adiabatic; 
short dashes: perpendicular-field adiabatic. Top panel: spine (on- 
axis); bottom panel: inner edge of shear layer. The adiabatic mod- 
els have been arbitrarily normalized to match the SSL model pro- 
files at the inner and outer boundaries. 



A second possibility which has frequently been dis- 
cussed in the literature is that the flaring point marks the on- 
set of turbulence, or the position at which Kelvin-Helmholtz 



1982 



199£ 



)ilities become non-linear (e.g. Baan 1980 




Bcgclman 


Bicknell 1984, 1986; 


De Young 1996 


1 




sen et al. 


[Rosen & Hardee 2000 


)• 



We will show elsewhere (Laing 



conservation-law analysis favours the hypothesis that the 
flaring point is associated with a stationary shock, primar- 
ily because it suggests that the jet is over-pressured at the 
beginning of the flaring region. This does not, of course, 
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Discontinuity 



Figure 21. Grey-scale images of the rms magnitudes of the mag- 
netic field components as fractions of the total field for the SSL 
model. Top: radial component (B 2 ) 1 / 2 /B; centre: toroidal com- 
ponent (B 2 } 1 / 2 /^; bottom: longitudinal component (B 2 } 1 / 2 /^?. 
The arrows marks the location at the edge of the jet where the 
three field components are roughly equal, as discussed in the text. 
The radial component is constrained to be zero for the spine and 
the s = streamline in the shear layer. The values for the inner 
region are poorly determined (Table [?]) . 



exclude the subsequent development of entrainment (and 
presumably turbulence), as we now discuss. 



5.3 Evidence for interaction with the surrounding 
medium 

It is generally accepted that jets in FRI radio galaxies de- 
celerate by picking up matter, but it is by no means clear 
whether the prin cipal source of additional material i s mass 
loss from stars (Phinney 1983; Komissarov 1994; Bow- 



if: 



man et al. 1 9961) or entrainment across the jet bou ndary 



( |Baan 198p| ; |Bicknell 1984 |l98fi| ; |De Young 1996| ): both 



Layer of 
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Figure 22. A sketch of one possible geometry for the inner re- 
gion and its transition to the flaring region, incorporating a slow 
boundary layer which does not persist at large distances from the 
nucleus. 



tre. These gradients are prima facie evidence for interaction 
between the flow and the external medium. There is no rea- 
son why mass input from stars should generate such gradi- 
ents (Bowman et al. 1996), although a pre-existing gradient 
might be preserved as a jet becomes mass-loaded. The form 
of the transverse velocity profile in our best-fitting models 
varies surprisingly littl e as the jet decelerates, but the er- 
ror analysis of Section 4.2 shows that the situation might 



are expected to be important. Our models require significant 
transverse velocity gradients, in the sense that the edge of 
the jet is travelling about 30% more slowly than the cen- 



be more complicated: a top-hat velocity profile at the inner 
boundary is consistent with the data, so the profile could 
still evolve significantly along the jet. The presence of large 
quantities of very slow material at the edges of the flaring 
and outer regions is firmly excluded, however (Table [?]). 

A second piece of evidence favouring deceleration by 
interaction with the external medium is the complex field 
structure in the flaring region, where we were forced to in- 
troduce a significant radial component, increasing towards 
the edge of the jet, in order to explain the low degree of lin- 
ear polarization. This radial field component would not be 
expected from simple passive evolution of a mixture of lon- 
gitudinal and toroidal field in the smooth velocity field we 
assume. The most natural way to generate such a radial field 
component is for the flow to have a disordered, turbulent 
character towards the jet edges such as might result from 
large-scale eddies. This is precisely the situation expected 
at the edge of the jet in the initial "in gestion" phase of the 
entrainment process (De Young 1996). The velocity field is 
then likely to have significant small-scale structure which is 
not included in our model, but our estimates of average bulk 
flow speed are unlikely to be seriously affected. Even if there 
is no dissipation or dynamo action in such a turbulent flow, 
there will be significant amplification of the magnetic field 
by shear, so the sim plest adiabatic models, which assume 
laminar flow (Section 5.4), will be inappropriate. 

Another way to distinguish stellar mass loading from 
entrainment across the jet boundary is to ask whether stel- 
lar processes can provide the mass input rate required to 
produce the observed deceleration. It is clear from the work 
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of Komissarov (1994| ) and Bowman et al. (1996) that a jet 
which is decelerated purely by stellar mass loading will tend 
to reaccelerate on large scales, where the stellar density be- 
comes low but the outward pressure gradient and buoyancy 
force are still appreciable. Our models require continuous 
deceleration in the outer region, favouring boundary-layer 
entrainment as the dominant mechanism there. We address 



this question via a conservation- law analysis in Laing & Bri 
die (2Cj02"|), where we conclude that entrainment dominates 
after the beginning of the flaring region. 

Little is known about the properties of turbulent rel- 
ativistic shear layers, or of the viscosity mechanisms likely 
to predominate in magnetized relativistic jets. We cannot 
therefore relate the deduced velocity profiles to the inter - 
nal physics of the jets. We note however that Baan (198C) 



computed steady-state models for viscous jets in constant- 
pressure atmospheres and estimated both the transverse ve- 
locity profiles and appearance of the jets (on the assumption 
that the emissivity is directly proportional to the viscous 
dissipation) for several forms of the viscosity. Baan's mod- 
els generally predicted extended low- velocity wings that do 
not match our derived profiles. He did however discuss cir- 
cumstances under which flat-topped velocity profiles such as 
those inferred here might arise, including that of an electron- 
positron jet. 



5.4 The emissivity profile and adiabatic models 

We have determined the variation of noB 1+a (proportional 
to the emissivity) in the rest frame of the emitting material. 
Separation of this variation into particle and field contribu- 
tions requires additional assumptions. The X-ray emission 
from the jets (Hardcastle et al. 2002) is most likely to be 
synchrotron, rather than inverse Compton radiation, so we 
cannot use it to decouple the particle and field components. 
We therefore postpone a discussio n of the variation of pres - 
sure and density along the jet to Laing & Bridle (2002), 



where we also consider X-ray observations of the surround- 
ing hot gas. 

A number of authors (Baum et al. 1997, ^eretti et al 



199£; Bondi et al. 2000) have recently re-opened the possi- 



bility that the jets in FR I radio gala xies are adiabatic in the 



sense first defined by Burch (1979), i.e. that the particles 



suffer only adiabatic energy losses, there are no dissipative 
processes causing particle acceleration or field amplification 
and the magnetic field is convected passively with the flow. 
We defer a full discussion of this question to a later paper, 
since our data and models are both substantially more com- 
plicated t han is allowed by t he analytical approaches in the 
literature ( Baum et al. 1997 ) . The simplest adiabatic models 



do not allow for any turbulent flow (Section 5.3) and there 



is independent evidence for p article acceleration in 3 C31's 
jets from X-ray observations (Hardcastle et al. 2002). Nev- 
ertheless 
points. 



we can make a number of preliminary qualitative 
We take the analytical formulae from 



Baum et al 



(1997). In the absence of velocity shear and in the quasi- 



one-dimensional approximation, the field components vary 
as: 

Hi oc x 

B t oc {xpvy 1 



B r oc {xfiVy 1 

where x is the jet radius. For a purely longitudinal field, this 
leads to a variation of the rest-frame emissivity: 
e oc (r^-^+BVs^-Cioa+iaya 

= (r/3)- 1 - 37 ^- 5 - 83 

and for a perpendicular field (Bi = 0): 

(5a+6)/3 -(7c+9)/3 



(17?)- 
(vpy 2 - 92 x 



X 
-4.28 



The inner region poses a severe problem for the simplest 
adiabatic models: we have no evidence for deceleration so, 
if the conical region is fully filled, we would expect a very 
rapid brightness decline away from the nucleus (oc z~ 4 2S on- 
axis even in the perpendicular-field case) compared with our 
estimates of oc z -1 ' 96 for the spine and oc z~ 1,33 for the shear 
layer. Even for the steepest emissivity fall-off allowed by our 
error analysis (Table [?]), the indices are grossly discrepant. 
We have already argued that much of the emission in the 
inner region may come from a surface layer (Section ^jj) and 
the assumption that the radiating material expands with 
constant opening angle may be invalid. 

For the flaring and outer regions, we have computed the 
emissivity variations for the parallel- and perpendicular-field 
cases using our model for the radius and velocity of the jet. 
The results are shown in Fig. |20| where we have normal- 
ized the adiabatic models to the observed emissivities at 
the beginnings of the flaring and outer regions. Two exam- 
ple streamlines are shown for the SSL model: on-axis in the 
spine and at the inner edge of the shear layer. The adiabatic 
models predict emissivities which fall far more rapidly than 
is observed in the flaring region: the deceleration is too little 
and too late to compensate for the expansion. In the outer 
region, by contrast, the perpendicular-field adiabatic model 
predicts emissivities fairly close to those observed. 

There are two other fundamental problems with the adi- 
abatic models. First, the field structure in the flaring region 
is not consistent with passive convection in a smooth, ax- 
isymmetric velocity field. Our assumed velocity field acts so 
as to shear an existing radial component, thereby amplifying 
the component along the flow. It cannot, therefore, create 
the region of approximately isotropic field at the edge of the 
flaring region starting with what is essentially a mixture of 
toroidal and longitudinal components. Second, the assumed 
velocity field cannot change the ratio of radial to toroidal 
field. It is clear from Fig. Ell that the radial component es- 
sentially disappears at some point after the flaring region. 

We conclude that simple adiabatic models could not 
describe the inner and flaring regions, even if more realistic 
field configurations and the effects of velocity shear were to 
be included, but that a model of this type may apply to the 
outer region, at least if the radial field component is mostly 
eliminated by the outer boundary. Further investigation of 
this set of problems is outside the scope of the present paper 
and will be presented elsewhere. 
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Figure 23. The best-fitting model of the 3C 31 jets viewed at various angles to the line of sight with a beam of 0.75 arcsec FWHM. 
Left panel: logarithmic contours with fixed sensitivity, i.e. with the same lowest contour in all plots. Right panel: logarithmic contours 
with fixed 750:1 dynamic range i.e. with the same percentages of the peak intensity in all plots. Both sets of plots cover ±27 arcsec from 
the nucleus and the angular scale is indicated by the bar at the top of the diagram. 



6 3C 31 AT OTHER ANGLES TO THE LINE 
OF SIGHT 

The best-fitting spine/shear-layer models require the jets to 
be at 52° to the line of sight. 

Fig. shows the appearance of these models if ob- 
served with the jet axis at other angles to the line of sight 
at a resolution of 0.75 arcsec. We have not modelled the core, 
but need to make a crude estimate of the dependence of its 
flux density on 9 in order to illustrate the effects of observ- 
ing with limited dynamic range. For these calculations the 
effective flow velocity in the core has been arbitrarily set at 
— 0.95 (r = 3.2) and its flux is assumed to scale with an- 



gle according to the predictions of a simple single-velocity 
model for a pair of oppositely directed, identical jets: 

S c oc [r(l - /3cos0)r (2+o,o) + [r(l +0cos9)]~ (2+ac) 

The spectral index of the core, a c = 0. The models are 
shown for 9=90° (jet axes in the plane of the sky) , then for 
9 decreasing in 20° steps to 30°. The final model is shown 
at 6=18° as this is close to the limiting case that our code 
can compute, wherein the line of sight lies inside the widest 
cone angle subtended by the jet outflow at the nucleus (in 
the flaring region). 

The left panels show how the jets would appear if they 
could be observed at these angles to the line of sight with 
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Figure 24. The best-fitting model of the 3C 31 jets viewed at various angles 9 to the line of sight with a beam of 0.25 arcsec FWHM. 
Left panel: logarithmic contours with fixed sensitivity, i.e. with the same lowest contour in all plots. Right panel: logarithmic contours 
with fixed 2048:1 dynamic range i.e. with the same percentages of the peak intensity in all plots. Both sets of plots cover ±10 arcsec from 
the nucleus and the angular scale is indicated by the bar at the top of the diagram. 



the same limiting sensitivity as in our VLA data for 3C31, 
using logarithmically spaced contours. Fig. |24| shows a simi- 
lar display for an observing resolution of 0.25 arcsec FWHM, 
emphasizing the changes in appearance of the inner jet and 
the start of the flaring region. 

The #=90° case is, of course, symmetrical with two iden- 
tical centre-brightened jets that lack w ell-defined intensity 
maxima at their ba ses. Images of 3C 44 9 (Feretti et al. 199£ ) 
and PKS 1333—33 (Killeen et al. 1986) show precisely these 
features. 



At 8=70°, we see the effects of moderate differences in 
the Doppler boost between the two sides: the base of the 
main jet appears brighter, and the counter-jet is both dim- 
mer an d less centrally-peaked t han the main jet (compare 
3C296; " 



Hardcastle et al. 1997). At 0=50° 



we see essen- 
tially the symmetries observed in 3C31's jets. Note that 
the counter-jet brightens on an absolute scale at lower in- 
clination angles 9 because each line of sight now intersects 
a longer absolute path length through both jets. At 0=30° , 
the bright base of the main jet is effectively contiguous with 
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Figure 25. The best-fitting model of the 3C 31 jets viewed at various angles 8 to the line of sight. Vectors whose lengths are proportional 
to the degree of polarization and whose directions are those of the apparent magnetic field are superimposed on selected total intensity 
contours. Left panels: 0.75arcsec FWHM (±27arcsec from the nucleus); right panels: 0.25arcsec FWHM (±10arcsec from the nucleus). 
The angular and polarization vector scales arc shown at the top and bottom of the figure, respectively. 
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the unresolved core, and at #=18° the characterizatio n of 
the structure as a "jet" according to the usual criteria (Bri- 



dle & P jerley 1984 ) would be questionable at this resolution 



ity of FRI jets (i.e. before bending and other intrinsic e ffects 
become dominant). The results of Laing et al. (1999) sug- 



By 0=18°, the apparent flux density of the core has also in- 
creased from 0.03 Jy at 0=90° to 1.6 Jy. This makes it un- 
likely that the wide-angle emission from the outer layers of 
the jet would be detected except in observations specifically 
designed for high dynamic range (or low angular resolution) . 

The right panels illustrate this by plotting the same 
five models with logarithmically spaced contours at fixed 
percentages of the apparent peak intensity, to a limit of 
1/750 of the peak. The contours are chosen to match the 
appearance of the left panel for the 0=90° case. These "con- 
stant dynamic range" displays probably correspond better 
to "standard" radio astronomical observations that have not 
been specially optimized to detect faint broad features in 
the presence of strong compact components. For 0=30° and 
0=18°, most of the emission detected outside the compact 
core comes from close to the spine of th e approaching jet. 
Images of BL Lac objec ts such as 3C 371 ( Pcscc et al. 200l"l 
Wrobel fe Lind 1990Q and Mkn501 QConway fc Wrobe] 
199 5]) s| how qualitatively similar features, although the ef- 
fects of projection exaggerate deviations from axisymmetry. 

Fig. ^B] shows the variation of polarization with angle 
at two different resolutions. The relative separation of the 
parallel-perpendicular apparent field transitions in the main 
and counter-jets from the nucleus is a strong function of in- 
clination. For the main jet, the transition point moves away 
from the nucleus into the flaring region as drops from 90 
to 45°, despite the opposite effect of projection on the po- 
sition of the flaring point. As decreases still further, the 
transition moves closer to the nucleus again. In contrast, 
the field transition in the counter-jet moves monotonically 
closer to the nucleus as decreases, with the parallel-field 
region being essentially invisible for < 30° . The longitudi- 
nal apparent field at the edges of the jets also becomes less 
prominent as increases, and would be difficult to detect for 
ss 90° in observations with limited sensitivity. Both effects 
are inevitable consequence of the toroidal/longitudinal field 
structure. We expect them to be general features of FRI 
sources, testable for complete samples even where detailed 
modelling is impossible because of intrinsic asymmetry or 
low signal-to-noise. There is little published data on field 
transition distances, but 3C 296 (Hardcastle et al. 1997) 
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and 0755+379 ( Bondi et al. 2000[ ) indeed have transition 
points further from the nucleus in their main jets. 

Figs |2^-^j| explicitly demonstrate the possibility of gen- 
crating a variety of apparent FRI jet structures with the 
same physical model by varying the orientation and of uni- 
fying FRI radio galaxies with some classes of "one-sided" 
objects. They also illustrate the need for careful considera- 
tion of observational selection effects when analysing statis- 
tical properties of unified models. In the presence of a range 
of flow velocities both along and across every jet, obser- 
vational selection through limited sensitivity and dynamic 
range translates into velocity selection within the jets. A key 
byproduct of our models may be a way to guide the statisti- 
cal interpretation of jet velocities in blazar -FR I unification 
models, as discussed by Laing et al. (1999). 



The analysi s of jet sidedness an d width ratios for a com- 
plete sample by Laing et al. (1999| ) shows that our models 
are likely to apply in detail to the inner parts of the major- 



gest that some model parameters vary systematically from 
source to source: in particular, the length of the inner re- 
gion and the characteristic scale of deceleration appear to 
increase with radio luminosity 



7 SUMMARY AND SUGGESTIONS FOR 
FURTHER WORK 

7.1 Method 

We have shown that an intrinsically symmetrical, decelerat- 
ing relativistic jet model containing simple prescriptions for 
the velocity field and emissivity with a locally random but 
anisotropic magnetic field, accounts for the major features 
revealed by deep VLA imaging of the straight segments of 
the jet and counter-jet in 3C31. The principal new features 
of our approach are: 

(i) the use of three-dimensional (but axisymmetric) pa- 
rameterized models of velocity, emissivity and field ordering; 

(ii) rigorous calculation of synchrotron emission, includ- 
ing both relativistic aberration and anisotropy in the rest 
frame; 

(iii) fitting to images with many independent data-points 
in linear polarization as well as total intensity using a robust 
optimization algorithm. 



7.2 Principal regions of the jets 

A major result of this modelling is that the three regions of 
the jet that were initially identified purely from the shape of 
the outer isophotes (Fig. are also regions with distinctly 
different internal variations of velocity and emissivity^]. 



7.2.1 The inner region (0 to 1.1 kpc) 

Our conclusions for this region are tentative because of the 
limited transverse resolution of the data. The region is char- 
acterized by: 

(i) low intrinsic emissivity; 

(ii) slow lateral expansion (a cone of intrinsic half-angle 
6.7°) and 

(iii) a significant component of emission arising in slow- 
moving material. 

The fitted central velocity is 0.8 - 0.9c. We have no evi- 
dence for deceleration in this region, but we cannot exclude 
the presence of higher-velocity (Doppler-hidden) flow com- 
ponents. Simple adiabatic models are grossly inconsistent 
with the emissivity profile. 



All distances in this section are measured in a plane containing 
the jet axis, corrected for projection using the angle to the line of 
sight for the best-fit SSL model. 
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7.2.2 Flaring region (1.1 to 3.5 kpc) 



7.3 Implications for unified models 



This region was defined initially by the more rapid spread- 
ing of its outer isophotes. Our modelling shows it to be a 
region in which several dramatic changes in the other jet 
characteristics occur together: 

(i) The jets decelerate rapidly to an on-axis velocity of 
0.55c after an initial slow decline from 0.77c. 

(ii) They maintain a transverse velocity profile in which 
the edge velocity drops to approximately 70% of the on-axis 
value. 

(iii) The intrinsic emissivity increases abruptly at the 
boundary with the inner region, then declines with distance 
from the nucleus, z, as z~ 3,1 in the shear layer and z~ 2 ' 5 in 
the spine. 

(iv) The emissivity at the edges of the jet drops to about 
20% of that on the jet axis. 

(v) The radial component of the magnetic field in the 
shear layer becomes significant, rising from zero at the spine 
boundary to 90% of the toroidal and longitudinal compo- 
nents at the outer edge of the layer, i.e. the field is essen- 
tially isotropic at the outer boundary of the shear layer in 
this region. 

(vi) The ratio of longitudinal to toroidal field strength de- 
creases slightly from about 1.1 to 0.8, independent of radius 
in the jet. 

The sudden increase in rest-frame emissivity at the flaring 
point suggests that there is a discontinuity in the flow, per- 
haps a stationary reconfinement shock system. The bright- 
ness and polarization structure in this region cannot be de- 
scribed by a simple adiabatic model. The transverse velocity 
profile and the growth of the radial field component strongly 
suggest that entrainment across the jet boundary becomes 
important. 



7.2.3 Outer region (3.5 to 12 kpc) 

In this region, the jets continue to expand on a cone of in- 
trinsic half-angle 13.1°. 

(i) The jets decelerate less rapidly, reaching an on-axis 
velocity of 0.26c by 10 kpc. 

(ii) The intrinsic emissivity in the shear layer declines 
more slowly (oc z~ 1A ) with distance from the nucleus. 

(iii) The transverse velocity and emissivity profiles re- 
main essentially unchanged from those in the flaring region. 

(iv) The ratio of radial to toroidal magnetic field strength 
decreases, becoming < 20% throughout the jet by 10 kpc. 

(v) The ratio of longitudinal to toroidal magnetic field 
in the shear layer continues to decrease, from 0.8 to 0.5 by 
10 kpc. 

Although the emissivity fall-off is much closer to that pre- 
dicted by a perpendicular-field, laminar adiabatic model, 
more work is needed to test this idea for realistic field and 
particle distributions. 

Beyond the end of the outer region, intrinsic environ- 
mental asymmetries begin to dominate, as evidenced by the 
large-scale bending of both jets. 



We have also calculated the change in appearance of our 
model brightness and polarization distributions as functions 
of orientation. These are in good qualitative agreement with 
observations of other well-observed jets and we therefore ex- 
pect the model (with some parameter variations) to apply to 
FRI jets in general. Figures |2^ and |2^ show that the inten- 
sity changes are considerably more complex than would be 
expected for single-velocity jets. They emphasize the need 
for high dynamic range and sensitivity to possible wide-angle 
jet structures when assessing whether observed jet proper- 
ties are consistent with unified models. We predict changes 
in polarization with orientation (Figure ^): these provide 
an independent test of unified models provided that our pro- 
posed field configuration is present in all FRI jets. 



7.4 Further work 

We now intend to model other resolvable bright jets in FR I 
radio galaxies to determine the extent to which their ob- 
served brightness and polarization properties resemble those 
of 3C 31. We expect to be able to infer their velocity, emissiv- 
ity and magnetic-field distributions, building on the broad 
success of the jet-deceleration model in accounting for the 
s tatistical asymme tries of the B2 sample of FRI sources 
(Laing et al. 1999). Other sources showing well-collimated 



inner jets and rapid flaring include NGC315 (Venturi et al 



1993), PKS 1333-33 (Killeen , Bicknell & Ekers 1986)arid 



3C449 dFeretti et al. 1999| ), and it seems likely that the 
regimes of collimation behaviour we have identified in 3C 31 
are common in FRI sources. We aim to study a sample of 
sources with a range of angles to the line of sight, if possible 
distributed isotropically, in order to test the results of Sec- 
tion [| We also plan to develop a more sophisticated error 
analysis in order to assess confidence levels with some degree 
of rigour. 

3C 31 has been cited as the archetypal FRI source, but 
is actually in the minority in having diffuse "tails" of emis- 
sion extending to large distances from the core rather than 
confined bridges an alogous to the lobes of FRII sources (De 



Ruiter et al. 1990). Significant differences in dynamics (es- 



pecially entrainment of the surrounding medium) might be 
expected between the two classes. We also expect that the 
deceleration process should depend on the jet power and the 

extern al environment. 

In Laing & Bridle (2002), we present a dynamical model 



for the jets in 3C31, based on the velocity field derived in 
the present paper, a description of the surrounding galactic 
atmosphere derived from Chandra and ROSAT observations 
( Hardcastle et a l. 2002j ) an d application of conservation 
laws following |Bickne!T*"(1994 ) . This approach should also be 
extensible to other sources. Our results favour entrainment 
across the boundary layer as the origin of the majority of the 
mass-loading of the jets in 3C31, but it will be important 
to explore this in other large-scale FRI radio galaxies. We 
should seek further evidence for the entrainment process, 
such as the reduced polarization near the boundaries of the 
flaring regions. 

Our ultimate goal is to replace the empirical descrip- 
tions of velocity, emissivity and field structure with realistic 
physical models. Although this is some way off, we have de- 
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veloped a self-consistent adiabatic model which can handle 
arbitrary field configurations and (laminar) velocity fields in 
a relativistic jet, with the aim of establishing whether any of 
the flow regions we have identified can be described in this 
way. 

If our interpretation of the emission from the inner re- 
gion of 3C 31's jets is correct, observations of the apparent 
brightness and motions of FRI jets on even smaller scales 
will not be sensitive to the properties of the underlying bulk 
flow, but only to those of its slowest-moving components, 
which may be essentially stochastic. Improved transverse 
resolution of the inner jets in such sources will be required 
to determine the origin and distribution of the slow-moving 
material, and the extent to which these innermost regions 
of FRI jets resemble the larger-scale jets in FRII s ources, 
e.g. those in 3C 353 (| wain, Bridle & Baum I998|). This 



will require greater sensitivity and longer baselines than are 
currently available with the VLA or MERLIN. 

Finally, a number of FRI sources (including 3C 31) have 
been detected at X-ray and/or optical wavelengths (e.g 
Hardcastle, Birkinshaw & Worrall 2001 ; Hardcastle et al 



2002; V /orrall, Birkinshaw fc Hardastlc 2001; |Sparks et al 



K 



raft et al. 2002; Perlman et al. 2001; Marshall et al 



200C 

2002 ) . ' ?he radiation is most plausibly produced by the syn- 
chrotron process over the entire observed frequency range, 
and the shape of the spectrum therefore carries information 
about particle acceleration and energy loss. It will be impor- 
tant to incorporate descriptions of these processes into our 
models. 
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